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ABSTRACT 

Some  efficient,  high  order  methods  are  discussed  for  approximating  the 
solution  of  an  initial  boundary  value  problem  for  a homogeneous  parabolic 
equation  with  time  dependent  coefficients.  The  methods  are  based  on  Galerkin- 
type  approximations  in  the  spacial  variables  and  single  step  methods  in  the 
time  variable.  The  equations  defining  the  time  stepping  procedure  are  solved 
only  approximately  however.  A preconditioned  iterative  technique  is  used  for 
this  purpose.  The  resulting  algorithm  is  shown  to  produce  optimal  order 
approximations  using  only  the  order  of  work  required  by  the  single  step  method 
applied  to  the  parabolic  problem  with  time  independent  coefficients. 
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Significance  and  Explanation 


Many  physical  situations  can  be  modelled  by  the  solutions  of  initial 
boundary  value  problems  for  partial  differential  equations.  Examples  of  such 


ituations  arise  in  the  theory  of  heat  conduction  and  other  diffusion  processes 


The  physical  parameters  involved  are  often  dependent  on  the  time  variable 


The  construction  of  accurate,  efficient  algorithms  for  the  approximate 


solution  of  these  parabolic  equations  is  studied  in  this  paper.  New,  effi 


cient  algorithms  are  proposed  and  completely  analyzed 
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EFFICIENT  HIGHER  ORDER  SINGLE  STEP  METHODS 


FOR  PARABOLIC  PROBLEMS:  PART  I 


James  H.  Bramble  and  Peter  H.  Sammon 


I . Introduction 

In  this  paper  we  study  efficient  ways  to  calculate  approximations  to  the  solution  of 
a parabolic  equation  that  are  of  third  or  fourth  order  in  time  and  of  high  order  in  space. 
The  approximations  are  generated  by  rational  function  based  schemes  (cf.  Nassif  and 
Descloux  [7]  or  Baker,  Bramble  and  Thom^e  [2]  if  the  operator  in  question  is  time  inde- 
pendent) but  these  schemes  are  modified  in  a manner  suggested  by  Douglas,  Dupont  and  Ewing 
[4]  in  their  work  on  the  Crank-Nicolson  method.  We  study  the  schemes  in  the  context  of  a 
linear  parabolic  equation  with  time  dependent  coefficients  in  this  part  of  the  work  and 
we  will  generalize  these  schemes  to  nonlinear  equations  in  Part  II  of  this  work. 

The  schemes  suggested  by  Nassif  and  Descloux  in  [7]  are  single  step  methods  based  on 

-x 

a certain  class  of  rational  function  approximations  to  the  function  e , x t 1R  and  a 
given  family  of  discrete  spacial  operators.  Nassif  and  Descloux  give  estimates  in  [7] 
that  show  that  the  resulting  approximations  make  errors  that  are  of  optimal  order.  How- 
ever the  schemes  are  not  really  suitable  for  practical  computation  since  each  step  of  the 
time-stepping  procedure  involves  the  solution  of  a new  linear  system  that  is  related  to 
the  family  of  spacial  operators. 

Douglas,  Dupont  and  Ewing  address  this  problem  in  [4]  and  suggest  the  remedy  of  using 
a preconditioned  iterative  technique  to  approximately  solve  the  linear  system.  This 
approach  only  requires  the  solution  of  linear  systems  involving  a fixed  discrete  spacial 
operator  if  sufficiently  many  iterations  are  done  at  each  time  step.  Conditions  are 
discussed  in  [4)  that  also  allow  one  to  iterate  a fixed  number  of  times  at  each  time  step 
(a  number  that  is  independent  of  the  discretization  parameters)  and  still  observe  the 
optimal  order  errors.  The  overall  work  required  by  this  strategy  is  of  the  magnitude  of 
the  work  required  by  the  usual  Crank-Nicolson  scheme  applied  to  a linear  problem  with  time 
independent  coefficients.  Thus  under  these  conditions,  they  obtain  a scheme  that  is  effi- 
cient as  well  as  effective. 

Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-75-C-0024 . This  material  is 
based  upon  work  supported  *•/  the  National  Science  Foundation  under  GrantsNo.  MCS76-07216  A02 
and  MCS78-09525. 


The  results  of  this  paper  are  similar  to  those  of  (41  with  regard  to  our  hiqher  order 
schemes  and  are  in  fact  in  some  ways  stronqer . In  particular,  because  the  schemes  which  we 
consider  are  inherently  more  dissipative  than  the  Crank-Nicolson  scheme , we  are  able  to  ob- 
tain the  best  results  unconditionally.  In  addition,  our  analysis  shows  that  the  special 
closeness  requirements  of  the  initial  values  to  the  "elliptic  projection"  demanded  in  |4] 
are  unnecessary  and  that  a more  natural  and  more  easily  computed  initial  projection  may  be 
used . 


An  alternate  approach  to  the  problem  of  finding  efficient  time  stepping  algorithms  ran 
be  found  in  work  by  Douglas  and  Dupont  in  [3] . They  analyze  two  efficient  schemes  for  para- 
bolic problems.  The  first  is  a method  which  is  of  first  order  in  time  (the  Laplace-modified 
procedure)  and  the  second  is  a three  level  method  of  second  order  in  time. 

We  now  introduce  the  parabolic  problem  and  some  convenient  notation.  Let  s)  c ik^  , 
d > 1,  be  a compact  domain  with  a sufficiently  smooth  boundary  3i!  and  an  outward  pointing 
unit  normal  n(x)  » (n^ (x) , . . . ,n^ (x) ) . Let  t > 0 . The  following  parabolic  problem  will 
be  studied  under  certain  smoothness  assumptions: 

r r 

-u  =*  L(t)u  3 - l D.(a.  (x,t)D,u)  + a (x,t)u  on  SI  x (0,i)  , 

i,j-l  13  3 ° 


(1.1)  / 
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o on  an  x (o,t) 


an 


Here  A = (a  (x,t)J . 


v on  n . 

is  a symmetric,  uniformly  positive  definite  family  of  matrices 


of  sufficiently  smooth  coefficients  on  n x a (x,t)  is  a nonnegative,  sufficiently 

o 

smooth  function  on  n x (0,r)  and  v(x)  is  a given  initial  data  function  on  i!  . If  the 
Neumann  boundary  conditions  are  under  consideration,  we  will  further  require  that  a (x,t) 
not  vanish  on  ft  x [0,t1  and  that  the  coefficients  Uv^}  have  the  following  special  form: 

a_(x,t)  = a(x,t)a^(x)  , 1 £ i,j  £ d , 

for  sufficiently  smooth  functions  a and  {a^}.  (This  extra  requirement  ensures  time 
independent  boundary  conditions).  We  will  refer  to  the  Dirichlet  boundary  conditions  as 


Pr,_  and  to  the  Neumann  conditions  as  BC. 
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o 2 

We  let  H denote  the  usual  L <a)~ based  Sobolev  space  with  norm  II  - II  ^ , 

negative  integer.  We  also  let  HX  denote  the  subspace  of  HX  that  consists  of  functions 

2 

that  vanish  (in  the  sense  of  trace)  on  3S1  . We  will  use  (•,*)  to  denote  the  usual  L (ID- 
inner  product  on  (1  . 

2 

The  operators  (L(t) torm  a fajnily  of  L (0) -self adjoint  elliptic  operators  on 
the  following  domain: 

C 2 1 

H n H if  we  have  BC„  . 
o D 

H2  n {w  c H2:nAVwl  =0}  if  we  have  BC  . 

I "Un  N 


Moreover,  the  form 

d 

D(t)  ( • , •)  - l (a..  D.(*),D,  (•))  + (a  (•),(•)) 

i,j-l  ij  j i 

is  (strongly)  coercive  over  H1  x H1  if  we  have  BC  or  H1  x H1  if  we  have  BC  . 

o o D N 

Thus  we  can  apply  the  standard  parabolic  equation  theory  to  (1.1)  (cf.  Friedman  [5]  or 
Lions  and  Magenes  (61 ) and  get  the  existence  and  uniqueness  of  solutions  u for  various 
classes  of  initial  data.  We  will  always  assume  that  v £ D and  further  smoothness  and 

li 

compatibility  conditions  will  be  added  later. 

We  let  T (t)  : L2(fl)  ■»  D denote  the  solution  operator  for  L(t);  that  isr  L(t)[T(t)fJ 
L 

» f for  all  f e L2U5).  We  note  that  <T<t>^<t<T  is  a smooth  family  of  bounded  operators 

from  H to  H n D^,  for  l > 0 and  that  {L(t))0<t<T  is  a smooth  family  of  bounded 

operators  from  Hi+2  n to  HX  , for  l >_  0 . In  fact,  (t)  = (-^-)^L(t),  j _»  0 , 

can  be  calculated  by  differentiating  the  coefficients  of  L(t)  with  respect  to  time  and  if 

we  let  T^  (t)  = (“■)  ^T  (t)  , we  have  that 
at 

T(1)(t)  * -T  (t)  L<X)(t)T(t)  . 

We  shall  use  the  symbol  C to  denote  a generic  positive  constant  throughout  this 
m2 

paper  and  we  will  define  £ (•)  = 0 if  m2  < . 
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1 1 . Hjs»«  lal  in  sol  et  l«at  ion  Opel  a ; >u  ■ 

We  will  assume  that  we  have  a finite  dimensional  subspacc  s.  ■ l * l .)  (associated 

h 

with  parameters  0 • h • l and  r ^ and  a sufficiently  smooth  family  of  sol f.-.d joint  , 

is'sitive  semidef  mite  operators  (T,  (til.  on  I "(si  that  have  range  in  S,  and  are 

h 0 <t  < i h 

positive  definite  on  . Ke  define  l^(t)  on  as  the  inverse  of  T^(  1 1 | ^ for 

each  0 t «_  i . Given  f • l.  (ill  , we  will  regard  T^ltlf  as  an  approximation  to  T(t!f. 
In  fact  we  will  require  that  the  following  he  ttue: 

(2.11  I (T*  (t)  - T^1' (tllfl  *_  C(j)hU2»fl  , ) ^ 0 . 0 1 1 1 r-2  i 


hi»rr 

<>:) 

* l.  (ill  . Finally, 

we  will  assume  that  there  is  a norm 

*'i 

on 

h 

t hat  Hat  inf  i«*s 

thi* 

f v'l  lowing  ; 

(2.21 

***3  < C*V»J 

* V.' 

(I^(t)v.vl  for 

0 ' t « i , y < s , 

- - h 

(2.  11 

|(l^(tl,i 

*2'  1 

10(1)1^1,1*^ 

for  0 . t ' t < S.  and  j ' 

— — 12  h — 

here 

<"  - 

^h  " 

Wo  not  e that  many  of  the  well  Known  oalerKin-tyi'e  met  hod  s satisfy  these  assumptions. 


For  a discussion,  see  Sammon  (S)  , 

We  let  P,(t)  T^(t)L(t)  iP  * , 0 ^ t ^ t , define  a family  of  "elliptic  proioc- 

2 r 

t ion"  operators.  We  also  let  PiL  (ill  • denote  the  orthogonal  L (ill-  projection  ont 

S.  . (Note  that  T,  - rr.pl.  If  W , t('*‘  i P.  for  seme  0 •.  ( ' r-2.  then 

h h h I,  — — 

(2.41  II w - Pwft  ' II w - P wll  < Ch*  * It wHI  , 0 < t ' t. 

- I - (*2  - - 

We  let  P,(’’(t)  - <-^-1  1P.  (t  1 for  j ' 0 . 

I ill  I 

Suppose  we  choose  v . H*  so  that  u(t)  . H*  and  llu(t)S  v ell vlt  for  0 «.  t < i . 

r - r — - 

Then  setting  W(tl  - P^(tlu(tl,  we  have  the  following: 

(2.M  lu(t)  - W (t  1 0 s c h'llvll  . 

— r 

w«'  wijih  to  how  will  t ho  t imo  ilirivAtivf*  of  w Approx iroat  o those*  of  \i  . To 

this  rml  wo  ntiuty  the  following: 


Proposition  (.'.I):  If  w 


. ( *2 


' P for  sv»i)  0 


r-2 , 


then 


Proof : 


lpjm  (t)wli  ^otmlh'*'  Hviy,  , 
Wo  see  that  (2.1)  implies  that 

«pT(m,(t)w«  < it  V t>^Vm-J)wii 

1 ~ j-o  1 h 


0 - t 


0 . 


m 


< II  l <™)T(1V"'-" 
j-0  ’ 


w II  * 0 h ( * ' II  wll 


II  (~ I™  TLvd  ♦ C h * * ' II  wll  - 0 h 


{ + 2 


II  wll 


t + 2 


An  easy  application  of  the  above  result  shows  that  if 


( ' r-2, 


( ♦ • ♦ «'w 

v < H is  suitably  chosen  then 


(2.h) 


llu(m)  (t)  - W,n'  (till 


C h<fJ  II v* 


t*2*2m  * 


an.l 


where  W(m  (t)  - Ip  particular,  this  shows  that  Klm'(t)  is  uniformly 

bounded  in  l.  (ft)  provided  that  u and  v are  suitably  bounded. 

At  times,  we  will  assume  that  the  following  condition  holds: 
iy  Hl^1’  (t)Th(s)ll  , #Th(s)l^j)  (DU  *C(j),  0 «■  s,  t v t,  ) M'  • 


Estimates  in  Sammon  1S1  or  Massif  and  nescloux  P)  show  that  this  condition  holds  foi 
various  Calerkin-type  met  Iks  Is  if  inverse  assumptions  are  valid.  The  following  result  is 
a consequence  of  Condition  : 

HU  (t)N<m'  (till  x V lli, * " (t ) II  for  m 2 I',  o x t . , . 
j-0 


-‘a- 


III.  Tim*  Piacr*t irat ion« 


r 


w*  now  consular  a method  of  oomput  mg  approximat  Ions  to  th*  aolul  ion  u(t>  v'f  U . 1 ' 

-x 

No  t>*g  1 n by  studying  rat  tonal  function  approximation*  to  th*  exponent  tal  * on 
IK*  . It  is  well  known  that  that*  at*  rational  function*  I ( x>  - 1 1*  and  0 at* 

relatively  prim*  polynomial*)  that  aatiafy  th*  following  condition*) 


(1) 

£(x)  ' 0 for  x ^ 

0 , g(0)  - 1 , 

().l)  y 

(ii) 

-!♦*  > \}*l(x)r(x)  ■>  l 

for  *>xtt»  >t  ' 0 , 

x ' 0 , 

(iii) 

. . -X | . *♦ l 

|r(x)  - * | ;C  « 

for  some  v M . 

x - 0 

. 

V 

wo  will 

ua*  r' 

a and  O'*  that  at*  no 

m>>t*  than  guadrat  lea 

In  our 

tat*r  veil  t>ut  for  now 

w*  assume  that 

( r*  < 

!»(*)  - i p. x Anvt  • ) 'l*x  wh*r* 

»\)"V 

l.  w*  l\av*  t h*  fol- 

<-0 

C-0 

l>>wing 

examples 

, where  r and  0 ar* 

of  d*gr**  two  or  l*»a: 

: 

(i) 

pi 

«0,  gt-I  with  >t«l, 

v-l 

(Hackwarda  Kttlei  > 

(it) 

*\m  ' 

1/3.  Pj-0.  g^-l/3 , g?-0 

with  \— 3 

(Ot ank-Nlcolaon) 

(lit)  Th*  family  patametrixed  by  l ' t/4,  \ * l/3t 

Pj-lJWt,  p,-ClJ-3l*l  'll,  g^-.M,  g,-l*  with  S - 0,  v - 
If  V-  * (Ul/.l)  . v-< 

Pj-p.-t’.  gj-‘.  g3-i/J  with  vt-i.i-.' 


( l v) 

(V) 

(vi) 


(t'alahan) 
(PadiO 
(Pad*) 
(radd) , 


Pj  — l/J,  p >*0,  g^-i.-M,  g.-l'b  with  >t  ' 0,  v * 1 
Pj-l'J.  p,-l  l-'.  gj-l /l.  with  >t  - 1.  v - 4 

w*  at*  particularly  interested  In  th*  oases  wh*i*  >t  - 0. 

W*  will  m>w  show  ?u%r  |'iv|vrty  ('.I'(iii)  can  l*ad  t>>  a two  point  Tayloi  expansion 
u«*d  by  Naaatf  ami  IVirltHia  in  I?). 

Proyosit  ton  (t.l)i  Snppo**  that  git'  t*  a amooth  function  on  (0,t  1 . Than  foi 
0 ^ t •.  f , • w*  hav*  t hat 

)»)),.,  V _ , ..  ).<)',.>,  . f*  I) 


()..') 


j"  g (-t)'g1'  (D  - j"  p (-t)  g ' (O'  ♦ / Klt.a'g1'* 
1-0  1 1-0  1 0 


(a'dn 


wh*r*  g'  (t)  * t“)'g(*'  and 


().  ») 


K it ,*) 


•I  . V • 


Proof : From  0.1)  (iii)  we  set'  that  \^(x) 


-x 


- P(x)  - 0(xV*S.  Now  take  m derivatives, 
where  0 ^ m ^ v,  of  each  side  of  this  equation  and  evaluate  them  at  x*0,  We  see  that 

m!rm  * l <™' (_1,m  1 3!  • 

m j-0  3 j 

This  gives  (3.2)  with  g(t)  “ tm,  0 < « m,  and  hence  (3.2)  if  g(t)  is  a polynomial 
of  degree  no  more  than  v. 

We  can  expand  a sufficiently  smooth  function  g(t)  in  a Taylor  series  of  degree  v 
and  apply  our  above  work.  This  shows  that  (7)  holds  with  the  kernel  given  and  completes 
the  proof. 

If  we  let  (u.  (t)  ) c s be  defined  by 

h 0<tv,  h 

(3.4)  u ♦ L.  (t)u  - 0 , u.  (0)  « v.  , 

h,t  h h h h 

where  v t S.  is  some  function  "close"  to  v ( for  instance,  v.  = Pv)  then  work  by 
h h h 

Sammon  (8)  shows  that 

(3.5)  llu(t)  - u.  ( t ) (I  - 0(hr)  , 0 - t < r, 

h — — 

under  certain  conditions.  Thus  if  we  could  approximate  the  solution  u^(t)  of  (3.4) 

with  a known  small  error,  we  could  use  (3.5)  to  show  that  our  approximation  is  actually 

close  to  u(t).  We  will  use  our  two  point  Taylor  polynomial  to  construct  an  approximation 

to  u^  . 
n 

Let  0 < k < 1 so  that  Mk  * r for  some  integer  MM.  We  will  study  a method  of 

P (x) 

approximating  u (k)  . Choose  a rational  function  r(x)  « -7  that  satisfies  (3.1)  (i) 
h Q(xl 

through  (3.1)  (iii)  and  where  the  degrees  of  P and  Q are  two  or  less.  (This  implies 
that  v <_  4.)  If  we  note  that 

uhU  ■ aK  ’ 'Vt,uh(t)  - 

uh2’  * ‘ {liM  * h!1,(t,,uh(t)' 

then  setting  g(t)  * u^(t),  0 ^ t k and  vising  (3.2)  gives  us  the  following: 

(3-f>)  (I  ♦ qjkL^tk)  ♦ q,k'(L‘(k)  - (k) ) )u^  (k) 

- (I  ♦ Pjkl^UD  ♦ r?k‘(L^(0)  - t^1’  (0))  'u  (0)  + 0(k''*1uh<V+1))  . 


-7- 


Thus  if  the  quantity  in  the  first  set  of  braces  is  invertible,  we  might  expect  the  fal- 


lowing to  be  "close"  to  u.  00: 

h 

2 2 (11-1  ^2  (1 ) 

V - il  ♦ q1kLh  ♦ q^k  (Uh  - 1^  )}  (k)U  * p^L^  ♦ p.,k  (L*  - 1.^  ) 1(0)  v,  . 

This  scheme  for  approximating  u^ (k I will  bo  the  basis  of  a scheme  for  approximating 
u^(Nk)  for  any  N 1 . It  will  be  seen  that  the  approximations  can  continue  to  be  de- 
fined using  operators  that  are  constructed  from  p , Q and  the  family  (l.^'ht)).  ^ . 

We  note  that  those  schemes  can  be  defined  if  the  degrees  of  r and  £ are  higher  than  two, 
as  was  done  in  [7) . 

The  solution  u^(t)  of  (3.4)  only  plays  a motivational  role  and  will  not  enter  in  any 

way  in  the  rest  of  this  work.  For  purposes  of  our  estimates,  we  will  need  some  function 

in  S that  is  uniformly  close  to  u(t)  in  the  sense  of  (3.5)  and  u.  (t)  would  be  a 
n h 

possible  candidate.  We  will  use  W(t)  for  this  purpose  however,  mainly  because  it  is 
easy  to  estimate  its  time  derivatives. 

For  0 ' n - M , let  t » nk  , L<3)  - L*3’  (t  ) ( j - 0)  , Tt3'  = T<3,(t  Hi  ' 0)  , 

~ — n n h n — n h n — 

2(1)  - 2(1) 

r - P(-kL  > , Q - ),  r » P - pjt  L , e - Q -q.k  L . Wo  now  cot  tie  the 

n n n n n*2n  *n  2 n 

question  of  Cn*s  invert ibility  on  . 

Propos i t ion  (3.2);  For  k sufficiently  small  we  have  that 

(3.7)  C,  (0  v*,v*)  v (o  v ,*')  < C_(0  * ,*)  , n 0,  **  t S,  . 

I n - vn  — 2 n — h 

Proof : This  follows  immediately  from  the  following  inequality: 

I 1*^2  ^n  ' — Ck(Qn*% ,v  ^ ’ 

We  will  assume  that  k is  sufficiently  small  for  the  rest  of  this  work.  Then  since 

is  invertible  due  to  £(x)*s  positivity  on  IK*  , the  above  Proposition  shows  that 

o is  invertible, 
n 

We  now  return  to  our  description  of  an  approximation  to  the  solution  of  (1.1).  Given 

0 0 
v(x)  , we  choose  a v « that  is  close  to  v (for  instance  V Pv  , although  we 

h 

n+  j 

will  describe  perhaps  better  possibilities  later)  and  recursively  define  V (0  - n - M' 
given  V°  , by  the  following  formula: 


P V , 
n 


0 < n < M . 


We  expect  V 


to  be  close  to  u.  (t  ) which  is  in  turn  close  to  u(t  ) (for  n > 0) 
n n n — 

and  we  will  derive  corresponding  estimates. 

2 

As  noted  before,  this  approximation  scheme  has  been  studied  in  an  L (ft)  setting  by 
Nassif  and  Descloux  in  [7] . We  now  see  how  computation  of  this  scheme  involves  solving 
a new  linear  problem  at  each  time  step  and  why  a more  efficient  variant  would  be 
desirable.  We  shall  later  study  the  variant  suggested  by  Douglas,  Dupont  and  Ewing  in 
(4] . However,  since  the  analysis  of  this  variant  requires  estimates  of  the  original 
scheme  (the  one  defined  by  (3.8))  in  new  norms,  we  shall  first  present  another  analysis 
of  this  scheme.  We  shall  also  define  a natural  choice  for  . 


IV.  Preliminary  Error  Estimates 


We  are  primarily  interested  in  how  close  vn  is  to  Wn  H W(t^)  (0  n M)  since 

we  already  know  (recall  (2.6))  how  close  wn  is  to  u°  = u(t^)  (0  £ n <_  M)  . As  noted 

2 

before,  these  estimates  are  already  known  in  the  L (SI) -norm  but  we  wish  to  study  them  in 


the  (possibly)  stronger  norms  given  by  (Cn (•),(•))  . This  will  allow  us  to  study  a 

1/2 

variant  of  the  scheme  where  the  (Q  (•)»(•))  norm  is  in  some  sense  a natural  norm  of 
the  problem.  We  note  that  most  of  our  work  will  go  on  in  in  the  next  two  sections. 

Letting  En  = vn  - Wn  for  0 £ n M,  we  see  that 

(4.1)  Q ^,En+1  = P . En  + (P  " P ^,)En  + (P  - P )e" 

n+1  n+1  n n+1  n n 

+ (Q  - Q .)En+1  - (Q  .,Wn+1  - P w")  . 
n+1  n+1  n+1  n 

This  will  be  an  important  error  equation. 

We  note  that  (4.1)  is  of  the  form  QW  = PV  + F where  W,V,F  e and  Q and  P 

are  selfadjoint  operators  on  that  satisfy  the  following: 

r 

(i)  (Q*,*>)  > 0 , 0 / <f  e Sh  , 

(4.2)  / (ii)  > 0 , 0 ^ £ Sh  , 

(iii)  ( (Q+P)?,ifi)  > HQf  ,<f)  , 0 / v?  e S , for  some  6 > 0 . 

V.  h — 

(Of  course  (4.2)  (i)  through  (4.2)  (iii)  follow  from  (3.1)  (i)  through  (3.1)  (iii).)  This 
situation  leads  to  the  following: 

Proposition  (4.1):  Let  QW  = PV  + F where  Q and  P are  selfadjoint  operators  that 

satisfy  (4.2).  Then  we  have  that 

(4.3)  (QW,W)  £ (QV,V)  - $((Q-P)V,V)  + 2(F,W)  . 

Proof : We  see  that 

(QW,W)  = (PV,W)  + (F,W)  = (PV,Q_1PV)  + (Q-1F,PV)  + (F,W) 

= (PV,Q~1PV)  + 2 (F ,W)  - (Q-1F,F) 
and  that  (PV,Q_1PV)  = (QV,v)  - ( (O-P) (I+Q_1p) V,V) . 

A simple  calculation  shows  that  I+Q  1p  is  selfadjoint  in  the  inner  product  defined 
by  ( (Q— P ) • , * ) and  hence  has  real  eigenvalues.  Thus  if  u is  an  eigenvalue  associated 
with  the  eigenfunction  , (4.2)  (iii)  shows  that  y(Q'ft^)  = (Q(I+0  ^P)v>,v5)  = ((Q+P)<^,^) 

>.  6 (QV  ,'f)  . 

-10- 


We  conclude  that  the  smallest  eigenvalue,  u , satisfies  u > 6 . Again  using  the  self- 

m m — 

adjointness,  we  see  that 

( (Q-P)  (I+Q"1P)V,V)  > Um((Q-P)V,V)  >_  6((Q-P)V,V) 
which  completes  the  proof. 

We  will  want  to  apply  Proposition  (4.1)  to  (4.1)  and  obtain  an  estimate  for 
(Qn+^,En+*,En+^) . In  anticipation  of  this,  we  prove  the  following  estimates. 

Lemma  (4.2)  ; Suppose  n^O,  0 _<  m ^ C and  we  have  Condition  if  Q(x)  is  quad- 

ratic . Then 


(4.4) 


l((Vm  ' VW 

*(2n+m  - VW 


< C k(Q  ,V>.)1/2(Q  ,*,)1/2 

— nil  n 2 2 


Proof : We  first  note  that  (3.1)  (ii)  implies  that  the  degree  of 

the  degree  of  Q(x)  and  (3.1) (ii)  and  (3.1) ( i ii ) imply  that 
constant  1 . Also  if  we  let  R(x)  = 1+x  if  q2=0  and  R(x)  = 
see  that  R 1(x)Q(x)  , R(x)Q  1 (x)  <_  C for  x ^ 0 . Thus  letting 


we  have  that 


P(x)  is  no  greater  than 

Q(x)  cannot  be  the 
2 

1+x+x  if  q2  j*  0 , we 

Rj  = R (kLj ) (0  _<  j <_  M)  , 


C (R  < {<}.*,*)  < C,(R. *>,*>),  £ S.  , j > 0 . 

l]  3 ~ ^ 1 h — 


Thus  it  suffices  to  prove  (4.4)  with  R^-inner  products  on  the  right  hand  side. 
We  have  the  following  estimates: 


(4.5) 

| ((L  -L  >*.  ,y>  ) 
n+m  n 1 2 

| < C k sup 

It  -si  <mk 
1 n l — 

KLh1>(s,W  1 

< CklkjH 

II*,  llT  ( Ck  (L  f ,' 
2i  — n i 

. )1/2(l  - )1/2  , 

l n 2 2 

(4.6) 

1 < 

n+m  n 1 2 

’"WW 

L .JPj)  I + | (L  *.  , (L  -L  )V>_)  I 

n+m  2 n 1 n+m  n 2 

< C k ( sup 

llL*1)  (s)T  II  ) ( 1-4-11  L T II ) II  L V,  II II  L *,ll 

h n n+m  n n 1 n 2 

It  -s | < m k 
1 n ' — 

note  that  Condition  B was  used  in  (4.6). 
n 

We  now  can  use  (4.5)  and  (4.6)  to  complete  the  proof. 
Lemma  (4.3)  : Suppose  n ^ 0 , Then 


(4.7) 


< C k Vrv’l)1/2'V2^2,1/2 


-ii- 


Proof : 


Wo  note  that  there  is  nothing  to  prove  unless  £(x)  is  quadrat ic  and  that  it 


suffices  to  prove  (4.7)  with  K inner  products  on  t lie  right  hand  side,  as  was  observed 
1 n 


in  the  proof  of  Lemma  (4.2).  We  have  that 

2 i (1)  . . \ i -.2. 


k^l  <LllV,  ,✓.,)  | < Ck“ll*,ll  IIntJI,  -•  ck'  (I.  )l/2(L  *.,*,)' 

1 n 12'—  1 l 2 I — nil  n 2 2 

which  suffices  to  show  (4.7). 

We  now  study  the  truncation  error  term  in  (4.1)  by  comparing  it  to  the  true  solution 
u (t ) of  (1.1)  . 


>1/2 


Proposition  (4.4):  Suppose  that  v < It 


max  (2  ( v+1 ) ,r+2)  , is  such  that  II u ( t) II 


r+2  - 


Cllvll  , llu('+1)  (till  < Cllvll  for  0 < t < t and  we  have  Condition  B,  if  £)(x)  is 

r+2  — 2 (v+1)  — — h 


quadratic.  Then  for  0 ■ n < M and  V < S.  we  have  that 

— h 


(4.8) 


l<Qn+1"  • Wn'*>l  (hr|/vll^7  + kUllvll,f>iil,) 


1/2 


2 (v+1) ' 


Proof : We  note  that  u^  = -I,u,  u^ 

Proposition  (i.I)  implies  that 

n+1  . n+1  , 2 n+1 


2(1)  ( 1 ) 

(O-L  )u  ” -L(I+T  )u  for  0 < t < i and  that 


,.  „ , n+l  , n+l  .2  n+1.  . n . n , 2 n . „ „ ,v  +1  „ 

(4.9)  II  (u  _llikut  +ll2k  utt  * ~ <u  ■plklltfp2k  Utt  ick  11  v- II.,  (v+i)  • 

j 


We  also  have  that  W 


pi° 


pui  + (Pj-piu-1 


(4.10) 

(4.11) 


j > 0 , 


L.W  = L,  (t.)T.  (t.)PL(t.)  u(t.)  = -PU; 

3 h 3 h 3 3 j t 

(iZ-l/^lwi  - -L.  (P+T^hu^ 

1 j 3 1 t 

= -L. (PTU^+P  T(1) (t.)u^)  + L,  (P  -P)u:j 
j I t I ] t j I t 

+ L.  (P  -P)T(1)  (t.)u^  + L.P(T<1)  (t,)-T(1)  )uj 
j I 3 t j j 3 t 


j > 0 


and  -1. ^ Pj  + T*l'(t^)uj|)  = P uj?t  . 

We  now  use  these  facts  to  see  that 

tn+i 

|(Cn+lWnM‘PnWn^)  I — c kV+1  ,,vH2(v+1)H^11  + c(  / Hw(1)  (s)-Pu{1)  (s)llds)IMI 

fn 

+ C(ll  (P  -p)u"fl|l  + ll(P  -P)u"ll  + ll(P  -P)T(U  (t  ,)u"4l|l  + ll(P  -P)T(1,(t  )un|l 
it  I t I n+1  t I n t 

+ ll(TU,(t  .,)-T<,1!)u"fl||  + II  (T<1>  (t  )-T(U)un||)  k2  (II L ,*ll  + IIWlI) 

n+l  n+l  t n n t n + 1 n 

v Ck  (hrllvllr+2  + kVKvll  > <IMI‘  * k2m,n^ll2)1/2 


which  qives  our  result. 


We  can  now  put  these  insults  toqethei  and  demonstrate  a bound  In  a Q ' -norm  for 
th«>  <Uff*t«irp  between  v"  and  W1'  . 

Theorem  ( 4 . S ) : Sup|K>se  that  v is  sufficiently  smooth  And  compat  ible  (the  hypotheses 

of  Proposition  (4.4)  would  suffice)  And  thAt  we  have  Condition  Hj  if  ^> ( x ) is  quadratic. 


Then  for  0 * N M , 

(4.12)  lla,/2(VN  * WN) II  < C(hr  4 kU) II vll  4 Cll0yj(V°  - P v)ll  . 

N — u 0 I 

Proof  Let  0 ^ n • M.  An  application  of  the  results  of  this  section  to  the  terms  on  the 
riqht  of  (4.1)  shows  that 

(On+1En+1,En+l)  < (14C  k ) (QnEn,EB)  t Ck(kV  4 hr)‘  II  vll  2 , 

which  qives  the  result. 

We  will  now  examine  possibilities  for  the  startinq  function  V°  . We  require  that 
V - PjV  be  bounded  by  C(k't  hr)  in  the  Q^2-norm  if  we  wish  a comparable  error  in  (4.12) 
We  can  always  let  V*  • P^v  but  we  note  that  the  approximation  scheme  defined  by  (1.8) 
never  requires  that  we  determine  Tj  applied  to  any  function.  In  applications,  this  would 
amount  to  a special,  expensive  calculation  required  only  at  the  beqinninq.  There  is  another 
choice  for  V°  that  involves  solvinq  a system  with  Qq  (or  even  . Since  such  systems 

have  to  lie  solved  anyway  to  take  the  first  step,  this  would  seem  to  be  a better  approach. 

We  have 

Proposition  (4.6):  Let  v i H*  n D , 1.(0) v f D and  define  V0'1  and  V°’2  in  S,  by 

i*  L h 

the  following: 

I OeV0'1  - p(v  * q kl.(0)v  4 q,k2(L2(0)  - L(U(0))v) 

(4.1J)  y ’ „ 1 2 „ 


0.1  22  l\) 

I - P(v  t q kI,(H)v  4 q k (1,(0)  - l.  (0) ) v) 

.11)  y " , 1 2 _ 


P(v  4 q1kt,(0)v  4 q^k  L (0)v)  - P Q(kL(0))v 


(4.14) 

"oJ/2(v0'2 

Pr<X'f  : 

We  observe  that 

C h II  vll 


(L0P]V  - I.  (0)v,y)  - ((I  - P )L(0)v,t,  *) 


(L^’p-v  - L(1)|°lv,fl  - <<T<1)<0)  - T*n)L(0)v,l.  v")  4 ( (P  -I)T(”  (0)L(0)v,LV)  . 
n I 0 0 1 0 


-l  1- 


This  completes  the  proof  of  the  first  part  of  (4.14)  and  the  second  p>art  follows  even  more 
simply  from  the  above  observations. 

This  result  completes  our  error  analysis  for  the  approximation  scheme  defined  by 
(3.8)  with  V*'  defined  by  either  equation  of  (4.13).  We  will  call  this  the  base  scheme 
in  the  sequel. 

We  note  that  if  V^’2  is  defined  by 
(4.15)  Q(W.0)V°'2  - PQ(kL(0))v 

where  Q(x)  is  a polynomial  that  satisfies  0 >c  sup  5(x)/Q(x)  <"»,  an  estimate  like 

CXx<® 

(4.14)  will  still  hold.  Thismodification  might  prove  useful  if  £>(kLQ)  is  a preconditioning 
operator  for  the  kind  of  linear  system  solving  techniques  we  will  study  in  the  next  section. 
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V.  A Variant  of  t he  Bast*  Scheme 

As  we  noted  before,  the  calculat  ion  of  tin*  base  scheme  involves  the  solution  of  a 
new  linear  problem  at  each  time  step.  We  wish  to  study  a variant  of  this  scheme  where  wo 
only  approximately  solve  the  lineal  system  at  each  t ime  step.  We  pro|>ose  t o use  an  itera- 
tive technique  for  this  purpose  which,  as  we  will  see,  can  he  provided  with  a qood  initial 
guess  for  the  true  solution. 

If  we  are  at  a point  in  our  calculations  where  we  have  several  accurate  approxima- 
tions to  the  function  u(t)  at  previous  time  steps,  it  can  be  seen  that  there  is  an 
ext  ta|x>lat  ion  of  these  values  that  yields  just  as  qood  an  approximat  ion  at  the  next  time 
step.  The  smoothness  of  u(t)  makes  this  possib*«.  This  extrapolation  could  he  used  as 
an  initial  guess  for  an  iterative  procedure.  But  this  observation  raises  a question. 

Since  even  the  exact  solution  of  the  system  which  we  are  approximately  solving  Is  no  closer 
to  u (in  the  sense  of  order)  than  the  extrapolated  guess,  why  iterate  at  all?  If  we  did 
no  iterations  and  used  this  procedure  as  an  algorithm  to  generate  further  approximations, 
errors  would  grow  and  the  approximations  would  deteriorate.  Such  an  algorithm  is  not 
stable.  Of  course,  the  Ivise  algorithm  (solve  the  system  exactly  and  forget  about  iterations) 
can  easily  be  shown  to  be  stable  although  we  will  not  formally  state  this  result.  Also,  it 
will  not  be  too  hard  to  see  that  if  we  make  an  error  in  approximately  solving  the  system 
that  is  of  the  order  of  the  local  truncation  error  and  that  is  in  some  sense  independent  of 
the  initial  guess,  then  the  algorithm  is  stable  and  gives  accurate  approximations.  For 
the  iterative  schemes  that  we  will  consider,  this  strategy  requires  a quantity  of  iterations 
that  is  on  the  order  of  the  logarithm  of  the  total  number  of  t ime  steps,  per  time  step. 
However,  there  is  a more  efficient  strategy  available  if  the  polynomials  P(x)  and  £(x) 
have  the  correct  proport ies . If  one  does  in  fact  give  a good  initial  guess  t o the  iterative 
scheme  and  then  iterates  only  a certain  number  of  times  per  time  step  (a  number  that  is 
independent  of  the  total  number  of  steps),  then  even  though  accuracy  is  not  improved,  the 
resulting  algorithm  is  stable  and  generates  accurate  approximat ions  for  u . This  phenom- 
enon was  first  observed  In  141  in  relation  to  t hr  Oratik-N  tool  son  scheme.  We  will  give 
arguments  in  this  section  that  show  that  these  results  hold  for  schemes  that  have  the 
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right  kind  of  dissipation;  that  is,  P and  Q are  such  that  6 > 0.  Similar  results 
can  be  proven  for  polynomial  pairs  that  are  just  stable  < iS  = 0)  but  the  condition 
k <_  C h* , for  some  constant  C , is  required.  This  condition  introduces  dissipation  and 
was  used  in  (4). 

We  begin  by  discussing  the  properties  of  a particular  type  of  preconditioned  iterative 
technique  for  solving  linear  systems.  We  will  assume  that  we  are  workinq  in  a finite  dimen- 
sional space  H with  an  inner  product  (•,*)  and  a norm  IMI  « (•.•)X'/‘  • Suppose  A 

H H H 

is  a positive  definite  selfadjoint  operator  on  H and  we  wish  to  find  an  approximation  to 
the  vector  x that  satisfies  Ax  « y,  where  y is  known.  We  will  also  assume  that  the 
situation  is  such  that  we  have  another  positive  definite  selfadjoint  operator  A ^ at  our 
disposal  for  which  A^Xz,z  r H is  easy  to  find  and  for  which  we  know  the  following  spectral 
estimate: 


(5.1) 


X0(AoZ'Z)H  - UZ'Z)H  i VAoZ’Z)H  ' 


z < H , 


where  0 < X ^<  \ are  known  constants.  Then  there  are  methods  which  when  given  an  initial 
guess  x<0)  for  the  solution  x , generate  a sequence  of  approximations  x*"',  n > 1 , to 
x that  have  the  following  properties: 

(i)  The  calculation  of  x^a+X'  given  {x^)™_g  only  requires  evaluating  A and 
Ao<  solving  systems  involving  A^  and  Hilbert  space  operations, 

(ii)  The  sequence  x^  •*  x as  a •*  “ geometrically  in  the  following  way.  There  is 


a decreasing  function  0 y(£)  <1  (0  < £ 1)  that  satisfies  y(l)  = 0 and  which  gives 

.1/2 


the  rate  of  convergence  of  the  iterative  scheme  in  the  A norm: 

X ° 

(5.2)  II  a!' ‘<x-x '"')!..  < Cjy"  ( )||A1/2(x-x<0))II„  , u > 0 . 

' ^ O H 


,.1/2,  (ci) . „ _ a r 0 \ „ 1/2 , ( 

II A (x-x  )ll  < C.y  1 t — ) II A (x-x 
o H — A ' 7 o 

A given  method  may  or  may  not  actually  use  the  spectrum  estimation  constants  \ and 
in  its  calculations.  We  also  note  that  (5.2)  implies  that  the  following  estimate  holds: 

(5.3)  llAX/2(x-X(a))llH  < Cx(l-y“(  ~ ))"V(~)||Ay2(x(a,-x<0,)||H.  a - 0 . 

The  preconditioned  conjugate  gradient  algorithm  fits  into  this  framework  with 

1/2.  1/2 

y " (1-C  )/l+f.  i see  (1).  We  also  note  that  another  example  of  such  an  algorithm  is 


given  by  the  following  descent  method.  Let  p - 0 and  qiven  x 
define  x^'fX^  by  the  following: 


(a) 


for  some  a ' 0, 


-lb- 


f 


, <*♦!>  **  U> 

Ax  ■ v* Y ♦ (A  -wA)x 

o o 

■y 

This  method  converges  for  certain  values  of  y and  if  we  choose  y ■ - ‘ then  we 

1 vo 

have  a method  that  satisfies  (i)  , (ii)  above  with  > • (1-O/l  + t . 

We  intend  to  use  an  iterative  scheme  with  the  properties  described  above  to  approxi- 
mately solve  the  system  ('.8)  which  defines  our  hase  scheme.  The  Hilbert  space  11  will  lie 
with  the  L*  (t!) -inner  product  and  the  above  discussions  outline  (nissible  error  results. 

We  will  keep  the  conjugate  gradient  alqoritltm  in  mind  since  it  offers  a govxl  convergence 
rate  and  it  does  not  require  the  values  of  and  \ in  its  calculations;  they  only 

enter  into  its  error  analysis. 

We  now  must  decide  what  to  use  as  a preconditioning  operator.  We  note  that  the 

contribution  of  the  term  in  0 , is  small,  so  we  can  ignore  this  term  when  con- 

n*  i n* l 

structing  a preconditioner.  We  now  discuss  two  possibil it ies : 

A.  as  a precondit  joning  operator. 

This  is  a good  choice  if  linear  systems  Involving  are  easy  to  solve.  For 

instance,  if  k'(x)  " ltq^x  is  linear  then  - I+q  kl  has  essentially  the  same  structure 

*> 

as  the  L operator  and  solving  this  type  of  problem  is  well  studied.  If  fix)  - (l*\xl* 
as  in  the  Calahan  method,  then  solving  systems  with  Q - (I*-Vkl.  )*  only  involves  solving 
two  successive  systems  with  the  (1  * lkl.^1  operator.  Thus  o is  also  a goixl  preconditioner 
for  this  method.  If  fix)  is  not  a perfect  square,  the  fourth  order  diagonal  Fad^  scheme 
being  a notable  example,  there  are  other  methods  for  solving  systems  involving  £ that  use 
complex  arithmetic.  Thus  using  f as  a preconditioning  operator  is  possible  for  these 
methods.  We  will  offer  an  alternative  in  H however. 

We  observe  that  the  following  result  contains  (h.l)  for  0.  as  a precondit ioning 

0 

operator . 

Proposition  (1.1):  Let  0 m,n  v_  M and  assume  Condition  H if  fix'  i s quadratic. 


Then  for  i S.  , 
h 


U*c|tn-tm*k|)  (cnv,v)  v (i^c ,*•)  v (i*c|tn-tm*k|i  (tf 
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Proof  i We  see  that 


‘V-*1  ’ ‘V**’  4 4 (lSn-^nW^’ 

•ml,  ust no  th*  t schniguca  of  proposition  (4.21  and  (4.11,  that 

1(0  -o  W,^l  ♦ ((CL-01*.vH  I c|t  -t  .k|(Q  *,*)  . 
m n mm  n m n 

This  giv*«  th«*  second  inequality  and  the  first  is  done  in  a similar  fashion  usinq  Proposi- 
tion (1.2). 

Thus  when  we  are  using  the  iteration  technique  at  the  n-th  time  step,  n 2.  1 > Wl’ 

can  always  take  1^  • 1^'  - (lac  * and  we  can  expact  an  error  reduction  by  a factor 

of  at  least  >''<(1.C  t ,1  *1  after  a iterations.  Note  that  this  implies  that  y < Ct  , 
n«  l — iv*  l 

ao  that  y » C k for  the  first  few  steps.  Also  if  we  are  given  an  < > 0 , there  is  an 

a 1/2 

a 2 * ( independent  of  nl  so  that  > »_  t t 

Now  we  consider  another  poss ibi 1 i ty  for  a preconditioning  operator  that  is  useful  if 
<}(x)  is  quadratic  but  not  a 1'erfeet  square.  Let  V ' 0 and  set  S ■ I ♦ VkL^  for 
0 <_  n < M. 

I*.  « (ltlkl.^1  or  8*  as  a preconditioning  operator 

We  first  note  that  it  is  easy  to  solve  systems  using  these  operators)  that  is,  we 
only  need  to  solve  (perhaps  successive)  systems  with  the  (Ielkl.  1 operator.  We  can  also 
prove  the  following  result  by  the  methods  used  in  the  proofs  of  Propositions  (1.21,  (4.21 
and  (4.  t)  : 

Proposition  (A.  21;  let  0 - m,n  <_  M and  i be  the  degree  of  • Suppose  that 

Condition  h.  is  satisfied  if  C<x)  is  quadratic.  Then  for  „•  , S,  , we  have  that 
h n 

(S.51  C.lsV,*!  < (0»-.vl  < C,(sV,/l  . 

In  — m - 2 n 

Thus  when  we  arc  using  the  iteration  technique  .it  the  n-th  time  step,  n ^ 1,  we 
can  Always  take  \^/V^  • 0^/C^  l an<1  obtain  an  error  reduction  by  %i  factor  of  at  least 

yrt  after  a iterations.  If  we  are  given  an  i > 0 * there  is  an  a <,  C log(l/»  t ^ *)  so 
i\  1/2  i/2 

that  y • i t , whet  e we  assume  that  » t < V . We  cannot  assume  that  a is  inde- 
— n n 

pendent  of  n this  time.  However  we  note  that  if  we  do  i * c|loq(t  ) j ♦ 0 iterations 
• n — n 

at  the  n-th  time  step  where  n ranges  between  l and  some  N , N'M  , in  the  course  of 
caiculat  ing  an  appro*  i mat  ion  to  u(t  )#  then  the  average  number  of  iterations  per  time  step  is 


-IH- 


1 V M V 1 , 

N l \ “ N ( \ \ M ’ - c • 

n=l  n=l 

Hence  the  average  number  of  iterations  per  time  step  is  independent  of  N , for  largo  N . 

We  now  gather  these  ideas.  We  will  assume  that  we  have  chosen  a preconditioning 
operator,  which  we  will  call  PQ  and  we  have  Condition  if  Q(x)  is  quadratic.  Thus 

we  can  assume  that 

(5.6)  C,  (PQV1*)  < (C  *,*)  < C,(PQ^,*>)  for  0<n<M,  *>£S  . 

i — n z n 

We  also  assume  that  we  have  an  iterative  linear  system  solving  process  IP  which  uses  this 

preconditioning  operator.  We  now  wish  to  use  IP  to  calculate  approximations  to  the 

solution  of  Q^x  * F > assuming  we  have  been  given  the  right  hand  side  F , an  initial 

guess  v for  y and  a tolerance  6 > 0 . We  will  assume  that  there  is  an  a = a (6  ) 
Ao  . * n n n n 

“n 

so  that  if  x is  the  an-t^  iterate  of  the  process  IP  , then 

(5.7)  llV/2<X-X(  "Nil  1 SnllPC1/2(x-X0)ll- 

Finally,  we  will  make  an  assumption  about  the  total  number  of  iterations  needed  to  achieve 

1/2 

certain  tolerances.  If  6 = e t for  1 < n < N , where  e > 0 , then  we  will  require 

x N n n - 

that  — ) a < C for  large  N < M ; that  is,  we  only  need  finitely  many  iterations  per 

N *■*_  n — — 

n=l 

time  step,  on  the  average,  to  achieve  these  tolerances. 


With  the  process  IP  at  hand,  we  now  formally  state  a variant  of  the  algorithm  stated 

in  (3.8).  First  of  all,  given  a v(x),  we  choose  a U„  e S.  that  is  close  to  v and  given 

u n 

a set  (8  }M  , of  positive  tolerances,  we  define  Un+*  in  terms  of  {U"'  }?  „ , 0 < n < M-l, 

n n=l  j=0  — 

in  the  following  way.  We  use  enough  iterations  of  the  process  IP  (which  uses  the  precondi- 
tioning operator  PQ)  to  generate  an  approximation  lin+1  to  the  (true)  solution  [in+1  of 
the  following  system: 

(5.8)  Q ,5n+1  = P U , 

xn+l  n n 

where  the  error  made  is  to  be  less  than  the  tolerance  g ,,  in  the  sense  of  (5.7) . We  use 

n+ 1 
n 

(5.9)  Z (U)  = y y .U3 

n+1  jio  Yn+1,D 

for  certain  coefficients  (y  , . „ , as  an  initial  guess.  (We  will  fix  values  for  these 

'n+1,3  3=0  * 

coefficients  later  when  it  will  be  clearer  what  they  need  to  be . Of  course,  letting 
Zn+^(U)  = Un  is  a possibility  and  in  general,  we  will  never  use  more  than  the  past  few 
values. ) 


If  wo  redefine  K1' 


Un-wn  for  n 


wo  note  that  wo  havo  the  following  important 


identity,  an  analoque  of  (4.1): 
(5.10) 


0 - p k"  + u-  -p  ,)Kn  ♦ (rn-pVn 

n+l  n+1  n n+l 


n+  1 


♦ (C  ,-Q  , )E 

n+l  n+l 


is  y“  - ► «"> 

*n+l  n 


0 < n < M-l  . 


We  now  analyze  the  error  made  by  this  kind  of  approximation  algorithm.  We  will  begin 
by  studying  a result  that  is  easy  to  obtain  but  is  not  the  best  possible  for  our  situation. 

for  1 n <_  v-1 

(if  v 2)  and  to  an  error  of  B = k for  v < n < M.  We  note  that  these  latter  toler- 


We  will  briefly  assume  that  we  solve  (5.8)  to  an  error  of  B = k 

n 


n — — 

ances  imply  that  for  our  types  of  processes  IP  , we  must  do  on  the  order  of  log  (M)  = 
log(i/k)  iterations  per  time  step  in  general.  One  might  expect  these  tolerances  to  lead 
to  good  error  estimates  (if  we  use  the  appropriate  initial  guesses)  and  we  will  show  th.t 
they  indeed  lead  to  qood  estimates.  However  we  will  later  show  that  we  can  get  the  same 
type  of  estimates  for  a modified  algorithm  that  only  requires  finitely  many  iterations 
per  time  step,  on  the  average. 

We  will  assume  that  v(x)  is  sufficiently  smooth  and  compatible  for  this  discussion. 
In  particular,  this  implies  that  we  can  take  U11  to  be  an  approximation  generated  by  IT 
to  either  V13'*  or  V^'“  (recall  that  these  were  defined  in  Proposition  (4.6))  with  an 
initial  guess  of  zero  and  an  error  tolerance  of  8^  = hr  . (However,  recall  (4.15).) 

Our  algorithm  is  of  course  not  well  defined  and  in  fact  will  not  obtain  the  accuracy 
claimed  unless  we  make  some  special  choices  for  the  starting  guesses  required  by  the  process 
IP  . To  be  able  to  do  this  for  the  various  schemes,  we  introduce  some  specific  examples  of 
the  operators  discussed  in  (5.9)  as  follows: 

* 0 for  0 < n < M , Z(1?  (U)  = Un  for  0 < n < M , 

— n+ 1 — 

■ 2Un-Un  * for  1 < n < M , 

- 3Un-3Un_1  + I)"-2  for  2 < n >.  M , 

* 4Un  - 6t'n  ^ + 4l)n  - for  3 n < M , 

* 5t'n  - lOP"*1  + lOtl"-2  - 5Un“3  + Un-4  for  4 <_  n * M . 
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2n+°l(l,) 

Zn+21(U) 

Z°!  (U) 
n+l 

z<4) (U) 
n+l 


-» ( 5) 
"n+l 


(U) 


We  can  now  state  this  not  quite  best  possible  algorithm  in  its  entirety: 

Algor  it  Ian  (11-.  Use  11'  with  the  preconditioning  operator  !'f  to 

(1)  generate  an  approximation  U('  to  either  V11’*  or  V1  using  zero  as  an  initial 

1 J* 

guess  ami  a tolerance  $ * -7  h ; 

(2)  generate  an  approximation  Un*'  to  the  (true)  solution  of  (5.8)  using  J (U)  =>  Un 


as  an  initial  guess  and  a tolerance  8 


k in  the  range  1 < n+l  y-1  : 


(1)  generate  an  approximation  un+*  to  the  (true)  solution  of  (5.8)  using  j (U)  as 
an  initial  guess  and  a tolerance  8n+j  “ k , in  the  range  v _<  ntl  <_  M . 


Again  we  note  that  since  we  are  using  a tolerance  of  8 


k for  v < n+l  <'  M , we  need 


on  the  order  of  log(M)  « log(i/k)  iterations  per  time  step,  in  general. 


We  use  the  techniques  of  Section  IV  to  study  this  process  via  Equation  (5.10). 


We  first  observe  the  following: 


(5.11)  |<f  (Un+1-u"+1)  ,En+1)  | <_  ClIO1^  En+1|lllP01/2  (Un+1-Un+l) 


En+1|lpn+llleiu<un+1-2i+l(l’,)l1 

— CPn+l'l®n+l  En+1"("0^<Hn+1^;(E))ll 

+ "^+i(wn+1-  Znil(W,,l,)- 


(5.12)  Ilf  :(En+1  - Z,  (E))ll  < C V 

xn*l  n+l  — L 


j*n-i+l 


Hcj/Vll. 

3 


where  ial  or  v depending  on  n . If  we  have  Condition  B and  a suitable  v(x>,  then 

(5.13)  Ilf  (Wn+1-z(1’  (W))ll  < dlwn+1-z(i!  (Will  + C kllL  , (WIUl-Z(i!  (W))ll 

n+l  n-fl  — n+l  n+l  n+l 

< Ck'  ( sup  IIW<j)(s)ll  + sup  II l-h ( s ) W ( 3 5 (s)ll) 

0<S't  Ovsvt  , 

n+l  n+l 

0<j<i  0 <_i  ^i  - 1 

< C kl  II  vll  _ , . . , 


whore  again  i*l  or  v depending  on  n . 


'hus,  we  can  show  the  following: 


n M 

Thoor . 3 ) : If  wo  goner at o a sequence  of  approximations  {V  ^ using  Algorithm  (1) 

where  v * II"  (u  * max (2 (v+1) #2 (r+1) ) ) is  sufficiently  smooth  and  compatible  and  wo 

assume  Condition  B,  , then  for  N > 0 , wo  have  that 
h — 


(5.14) 


llCy2(UN-WN) 

N 


C (hr+kv)B  vll 


Proof:  Lot  n •>  0 . Wi  have  via  (5.10)  and  our  estimates  that 


(5.15) 


(Q  .E  , E ) • (UCk)(Q  E ,E  ) + C k II  vIP  (h  *k  > ' 
n*l  — n L 


+ C k 


n-1 

l 

i-max  (n-v* 1 ,0) 


(Q-E^.E3) 

0 


1 /2  n v r 

Also  if  0 ^ n <.  v-1  then  II Q E II  C(k  +h  ).  These  inequalities  and  (5.15)  give  our 

result . 

Thus  wo  have  optimal  order  errors  for  Algorithm  (1). 

As  a preliminary  to  analyzing  an  algor  i wbir.  that  requires  only  finitely  many  itera- 
tions per  time  step  on  the  average,  we  prove  some  results  for  the  following  situation. 


Suppose  that  a sot  of  approximations  {lP  )n  . , to  the  functions  (WJ}"  . _ are 

i=n-i+ 1 i=n-i+l 

given,  where  1 i <_  v+1  and  n 2.  i“l»  Use  the  process  IP  (with  the  preconditioning 

operator  Q)  to  generate  an  approximat ion  un+*  to  the  (true)  solution  pn+^  of  (5.8) 

using  Z**!(U)  as  an  initial  guess  and  a tolerance  8 , N 0 . 

n+1  n+1 

By  the  analysis  done  so  far,  we  already  know  the  following. 

proposition  (5.4):  Suppose  we  have  Condition  B and  v(x)  is  sufficiently  smooth  and 

compatible.  Then  for  any  r N 0 , we  have  that 

(5.16)  (?  En+l,F.n+1)  < ( 1+C  k ) (Q  En,En)  - 4((C  ~P  Ik'V") 

ntl  — xn  2 ln  n 

v r ^ n ?i-2  2 

+ C k ( (k  + h ) + 18  ,k  )llv!P 

n+1  h 


V ||C1/2(Ej-Fj'1)H2 
k . . . vn 

;i=n-i  + 2 


The  above  result  motivates  the  investigation  of  the  equation  stated  in  tho  following 
result . 


Proposition  (5.5):  Let  gw  = PV+F  where  g and  P are  selfadjoint  operators  on  . 

Then  we  have  that 

(5.17)  ( (g+P) (W-V) ,W-V)  + ((g-P)W,W)  = ((g-P)V,V)  + 2(F,W-V>  . 

Suppose  further  that  g and  P satisfy  (4.2)  (iii) . Then 

MQ(W-V)  .W-V)  + ((g-P)W,W)  <_  ((g-P)V.V)  + 2 (F  ,W-V)  . 

We  now  apply  Proposition  (5.5)  to  (5.10).  If  we  can  enforce  a certain  important 

condition,  namely  that  6 > 0 for  our  polynomials  P(x)  and  Q(x)  (recall  (3.1)),  we 

can  obtain  an  estimate  that  will  allow  us  to  analyze  the  last  term  in  (5.16) . 

Proposition  S.6)  : Suppose  we  have  Condition  , v e H'1  (u  = max (2r+2 , 2v+2) ) is 

sufficiently  smooth  and  compatible  and  5 > 0 . Then  there  are  constants  6*  > 0 and 

C.  > 0 so  that  if  6 , < 8*  , we  have  that 

o n+1  — 

(5-1**,  ♦ <<Qn+1-Pn+1)En+1,En+1) 

< (1+C  k)  ( (Q  -P  )En,En)  + Ck2(g  ,En+1,En+1)  + Ck2(gEn,En) 

— n n xn+l  n 


+ Ck2  ((h'tkV  + S2  k2l_2  }||  vll  * + C6^  . I llgA/^(E^-E^~X)ll^. 
n+1  0+1  j=n-i+2  3 

Proof : Our  usual  techniques,  with  different  uses  of  the  aritlimetic-geometric  mean  in- 

equality, yield  (5.18)  with  (1+C  k ) replaced  by  1 and  with  the  following  extra  term  on 
the  right  hand  side: 

(5.19)  f(((Q  ,-P  - (g  -P  ))E  ,E  , |. 

xn+l  n+1  xn  n n n ' 

2 

Note  that  g(x)-P(x)  = x + 0(x  ) by  the  accuracy  condition  (3.1)  (iii)  and 

2 

Q(x)-P(x)  > 0 for  x > 0 . If  we  redefine  R(x)  = x + |q_,-p_Jx‘  then  C^R(x)  <_ 

Q(x)-P(x)  <_  C^R(x)  for  0 < x < «°  . Thus,  under  Condition  , the  techniques  of  Section 

IV  show  that 

(5.20)  | ( ( {Q  ^,-P  x.)  - (g  -P  ))En,En)|  < Ck2  (L  En,En)  + C k3  |q  -p  lllb  En|l 2 

'n+1  n+1  n n ' — n f r2l  n 

< C k (R(kL  )En,En)  < Ck((Q  -P  )En,En)  . 

— n — xn  n 

We  could  now  combine  (5.16)  and  (5.18)  with  suitable  choices  for  the  parameter  i . 

We  would  then  essentially  find  that  llg3'/“EN|l  is  bounded  by  terms  that  are  0(hr+kv), 
terms  that  measure  the  initial  error  in  the  g3/"-norm  and  terms  that  measure  the  initial 
error  in  the  k (Q^-P  ) 3^‘-norm . The  projection  we  have  chosen  for  the  initial  data 


i „j-l. 


II 
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(as  given  in  Proposition  (4.6))  is  defined  so  that  it  is  computable  by  the  IP  process 

1/2 

and  so  it  leads  to  an  initial  error  that  is  good  in  the  -norm.  Unfortunately  it 

-1/2  1/2 

does  not  necessarily  lead  to  one  that  is  good  in  the  k (Q^-P^)  -norm.  We  could 
now  proceed  in  two  ways.  One  approach  is  to  let  be  P^.v  or  look  for  another  special 

approximation  which  is  good  in  all  the  required  norms.  But  since  the  process  IP  would 
probably  not  be  useful  in  generating  such  an  approximation  (the  spectrum  of  does  not 

bear  the  correct  relationship  to  the  spectrum  of  L^(0)),  a special  process  would  be  needed 
to  generate  only  U°.  Since  we  would  prefer  to  avoid  this  situation,  we  are  led  to  using 
(5.18)  in  some  other  way.  After  all,  it  was  only  the  direct  use  of  (5.18)  that  gave  this 
apparent  problem. 

We  have  the  following  result  which  combines  (5.16)  and  a variant  of  (5.18);  the  latter 
uses  multiplication  by  the  time  variable  to  avoid  potential  problems  at  time  zero. 
Proposition  (5.7):  Suppose  that  Condition  B holds,  v c HP  (u=max  (2r+2 , 2v+2) ) is  suf- 


ficiently smooth  and  compatible  and  6 > 0 . Then  for  each  e > 0 there  is  a 6**  > 0 

1/2 

that  if  £ min(e**,tn+1> , we  have  the  following  where  C is  independent  of  e ; 

(Q_.  1 En+1  ,En+1 ) + C— -^7—  II Qx/J  (En+1-En)  II  2 


+ c[  V((Vrpn+i)En+1'En+1>  -TT 


( (Q  -p  )En,En)] 
n n 


< (l+C  k)  (QnEn,En)  - (6/3)  ( (Q  -P  )En,En) 

n n 

+ Ck(h2r+k2v  + e2+1k2l"2)||vl|2 

n+i  u 


+ e_lTL  I HQ1/2(Ej-Ej‘1l||2  . 
j=n-i+2  J 


Proof : It  is  a rather  straightforward  computation  using  (5.16)  and  (5.18)  to  obtain  (5.21) 

with  6/3  replaced  by  5/2  and  a term 

Cc  ((Q  -P  )En,En) 
inn 

on  the  right  hand  side  where  > 0 is  arbitrary.  This  gives  the  result. 

We  can  now  define  and  state  results  for  our  final  algorithm: 

Algorithm  (2 ) : Use  IP  with  the  preconditioning  operator  to 

(1)  generate  an  approximation  U to  either  V ' or  using  zero  as  an  initial 

guess  and  a tolerance  ~ hr  . 


(2)  geneiate  an  approximation  un  * to  the  (true)  solution  of  (5.8)  w ing  *p' 

as  an  initial  guess  an«l  a tolerance  t<  . min(k  n,t;)  in  the  tango  1 nr  l 

nr  1 


where  8 * 0 if?  small. 

(t)  generate  an  approximation  l’M  * to  the  (true)  solution  of  (5.8)  using  (b) 

as  an  initial  guess  and  a tolerance’  8 , min(8,t'  ')  in  the  range  vr  1 nr  l M, 

nr  l nr  1 

where  8 x 0 is  small. 

We  note  that  ♦ *.»•  important  difference  between  this  algorithm  and  Algorithm  (1)  i- 
that  by  out  assumptions  on  the  ptoceas  IT  , we  only  have  to  iterate  a fixed  numbei  of 
time,  .it  each  time  -tep,  on  the  average.  Thus  we  are  demanding  fewei  iterations,  but  as 
we  will  soon  see,  we  get  the  same  convergence  rates. 

Theorem  ('■  8)  ; Suppose  that  Condition  1^  holds,  v * h'  (n=max  (rt  2 f2\>r2) ) is  suf- 
ficiently smooth  and  rom|vit  ible  and  ^ N 0 . Then  there  is  a (computable)  8 N 0 so  that 

n M 

the  approx imat ions  |P  ' goner.it  -i  * Algorithm  (2)  satisfy 

n -0 

IK''  ' (wf'-l>N)ll  ' c(kv+h*  lllvll  , N • o , 

N \\ 

N N 1 /•'  v*  * 

II W -H  II  « C t (k  th  ) II  v II  , N ' 0 , 

1 — N U 

llu(t  )-UN||  0(kVthr)llvll  , N " 0 . 

N — U 

N N 

Note  that  we  get  a superconvergence  result  for  W -U  in  the  II  • II  -norm  fot  times 
bounded  away  ft*™  zero. 

Proof : The  proof  is  almost  immediate  from  (5.21).  We  first  see  that 

(0  1 ,K,U  5 ♦ \\<3l(htnn-V")r  . C((Q  . -P  ll)Rntl.Fntl) 

nr  1 tit  1 tit  l nr  l 

(5.22) 

- C(h‘r4k-’V)llvl|2 
II 

fot  1 titl  v,  where  we  note  the  important  role  of  the  t variable  in  the  case  ti-0. 

-C  t 

If  N ___  vrl,  we  set  i v-t  l in  (5.;>l)f  multiply  the  inequality  by  e n<  fot  a 


suitable  c and  sum  over  t he  range  v « n N-l  . Then  by  choosing 


■ 0 suf 1 t . t ent  1 v 


small  and  using  (5.22)  for  0 n v,  we  obtain  the  first  restilt  . The  third  result  now 
easily  follow?*. 

The  remaining  usult  is  obtained  by  noting  that  if  c S , o ti  M and 

h — 

K(x)  x t |g.,-p.,  lx'  tlt.'ti 


Hell  c(t,  c,c)  - (R(M.  K\e)  • ' i (e  P )e,e). 

l n k ti  k n n 


I 


I 


V I Computat  ional  Considerat  ion?; 

Wo  conclude  this  paper  with  a fow  remarks  concerning  the  cevnput  at  iona  1 aspect-,  «-t 

Algorithm  ( 2 ) for  quadratic  p(x)  and  o(x)  . Suppose  that  the  function;.  U*  form 

til 

a basis  for  S.  . Moreover,  suppose  that  these  functions  have  tn»en  chosen  so  t li.it  1 meat 

systems  on  ttO  involving  the  matrices 

A - ((1,  v’.  , , A(U  - [ C I . ‘ 1 ’ v- . ,.'  > f!  . , , 0 m M . 

m m l 3 i #3*1  m m i j i, )*1 

g “ Kv«.  ,*.nJ  , . 

i j i,j-i 

have  acceptable  computational  properties.  Wo  now  wish  to  examine'  Algorithm  (.’)  in  this 


We  begin  this  discussion  by  making  some  definitions  and  ident i f icat ions . u*t 


0 < m < M.  If 


V * v? . c S 
L l l h 


G i"  f **j)  1 j1Bll  . 


)'  [G_1A  £1  .V  . 

j-1  1 


Thus  if 


5 ■ ZV) 


satisfies  O U = >f  then 


(I  + q.kG  ^A  t q k2(G  *A  G 'a  - G lA*"))r  - £ 
1 m *2  mm  m -* 


or,  equivalently, 

2 - 1 “>  ( i ) 

(6-3)  V (G  + V\,  f ^2k  V Am  - VA  >i 

' Ci  - • 

We  note  that  (6.3)  involves  a symmetric,  positive  definite  matrix  H . We  now  lot 

m 

be  one  of  the  polynomials  in  kl.^  discussed  in  Section  V.  bet  *q(x)  • 1 ♦ q^x 
„ 2 

f q,x*  define  its  coefficients  and  let 
2 

‘V  - G f qjkA0  ♦ q2k-’A0G'1A0  , 

which  is  also  a symmetric,  posit. ivc  definite  matrix  on  IK*  . Finally,  lowriting  ( r> . u) 
gives  us  the  following: 

(6.4)  < * ^il»H  * ^ ^ * C , * * ^ • n • IK* 

where  < •,•)  is  the  usual  inner  product  on  IK 

i i 

Thus  we  set'  that  S.  with  the  I. ’(*.’) -inner  product  and  11.'  with  the  mnet 

h 

product  are  unitarily  equivalent  under  the  identification  of  with  the  i-th  canonical 

j 

basis  function  in  IK  . This  implies  that  any  process  IF  on  S,  that  solves  sy  t cnv. 


1 


involving  v uuuj  * y.'  a pt  e«  "ini  1 t iouor  i torwally  * ••juivalent  t ■ {t-  • 

m 

on  U\'  th.it  solve.  systems  involving  « • usi.ihj  * 1 l<  .«  « 1 i • • .»nd  1 1 i ■ *nei 

m 

a I'Hhv',*;  will  u nally  bo  qivon  on  U<%l  in  |i  .ot  i,v,  thi  identili.  it.  si  ii,  I im  tin 

cot  i esix»nd i nij  proves  lb  on  s. 

h 

How  *1*'  wi  oomput e t h»»  I'oi't  f i * ' 1 1 *nt  ■ ■ {,  ' tot  th»  Im  ; expansion  ! t 

1 l i .1  , ' n M 

(un  t ion  > H,M ) , ,,  i } *V  1 defined  l»v  Abfontlun  The  linear  v tem  used  t •; 

0 _n  _M  ii 

putt*  U1,  lui*.  v*  l*f  a * .1  riqht  han«i  ide  when  t l.*  (:■’).  It  *-  ^ \ then 

i i 


(»*.  M 


Thu-;  '•  ' ‘.in  l*’  ,a  1 on  l at « l by  taking  mini  |M»iiurt  with  t and  an  bo  t.»und  I \ 


solving  a involving 


!’h«  1 iiioat  y.tom  u « l t,.  compnta 


n«  l 


t . >i  ■.  «m»  n 


has  l’  l)!  V-  v a ■»  right  hamt  ide.  It  we  know  ' (the  eoefti*  lent  |,<i  • M)  then 


(C  * P tkA  t rk‘(A  li  \ -A(1,))  n , 
I n . n n n 


(6.6)  r, 

no  that  ..  . m be  calculated  by  solving  one  ystem  involving  o.  and  > an  b«  ai.  n 


lated  by  solving  two.  Thus  we  t an  compute  the  t ight  hand  ide-.  of  the  vanoii  -\  tern  . 
It  we  know  { 1 ^ {)  • initial  guess  lot  the  iteration  pnnedure  i e»  t.  ibulo. 


Then  to  find 


n»  l 


via  tin'  process*  lb  , we  may  have  to  evaluat* 
-1  I 


l 


H , Atul  G ' * \\ 
!!♦  I 


solve  systems  involving  ‘ B and  do  various  Hilbert  • pa*  ••  'pci.it  ion  in  11-  with 


t lie  <i>,*  ) innei  ptodurt  . bvaluating  v,  to 


n»  l 


to  i l mi *| lit  1 oi  wai  \ . 01  <w 


cvei  , nv'te  that  ‘ a 1 in  l .it  i rnj 


\f On, -i n>  <*B+1vV  ’h ■ n;’ 

only  m«'. nr.  olvww  on.  sy.lm  inv.ilvin.i  . It  wo  hnvo  t..  I v<  t >,.•  •.  I om  1 1 l> 

it  y. 

In  IK  , we  observe  that  thi:-  is  Mjuivalent  to  olviug  B (Thl  b.  ctr.  n 

veil  lent  in  many  nituat  ions  as  b may  bo  easier  to  obtain  than  rhi-  wa  t ’ . i , 

in  («*.*>)  and  (<■*.«>)  . ) Wo  mu  .t  now  solv»  \ tom.  that  involve  i«.  11  .lx' 

for  aunt*  \ o for  instance,  then 

’V  (. ; • \kA  >o* 1 \ kA  ) 

O 0 

Thus  solving  ' Hn  U in  thi:.  c.t-.r  moan  solving  tw«'  vtom.  that  im  l\,  tin  mu  .1  i \ 

(»!♦  iKA^)  and  t'va  l n.»  t 1 n»j  » • orwc  . if  o ( x ) 1 not  a port  e»  t gum  e , • t * « 1 1 1 . 1 

could  be  used  to  efficiently  solve  system;  involving  B 

The  partieul.it  ehoieo  of  the  iterative  pi  00  will  deteimin*  u!ii>  1 >t  t ’ - 1’  -v. 

in  the  implement  at  1 on  .-r  Alaoiitlun 

-2«- 


oon  ; idoi  at  i.'ti  i . 1 rlov.ml 
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